APPLICATION NO. 09/826,118 

TITLE OF INVENTION: Wavelet Multi-Resolution Waveforms 
INVENTOR: Urbain A. von der Embse 



Copy of reference: 

Eigenvalue algorithm by Vaidyanathan and 
Nguyen entitled "Eigenf ilters : A new 
approach to least-squares FIR filter design 
and applications including Nyquist filters' 7 



l 



IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS, VOL. CAS-34, NO. 1, JANUARY 1987 



11 



Eigenfilters: A New Approach to Least-Squares 
FIR Filter Design and Applications 
Including Nyquist Filters 

P. P. VAIDYANATHAN, member, ieee, and TRUONG Q. NGUYEN, student member, ieee 



Abstract— X new method of designing linear-phase FIR fitters is pro- 
posed by nrfnfmnring a quadratic measure of the error in the passband and 
stopband. The method is based on the computation of an eigenvector of an 
aufupriate red, symmetric, and positfte-deflirite matrix. The proposed 
design procedure is general enough to incorporate both time- and 
frequetxy-domain constraints. For example, Nyquist filters can be easily 
rfrejgnerf using this mif™** The desi&i time for the new method is 
comparable to that of Remez exchange techniques. The passband and 
stopband errors in the frequency domain can be made eouhipple by an 
iterative process, winch involves feeding back the approxiniation error at 
each iteration. Several numerical desimi examples and comparisons to 
existing methods are presented, which demonstrate the usefulness of the 



N-l is even [7].) The desired response is 

U 0 < <o < 



I. Introduction 

IT IS WELL KNOWN that most linear-phase finite 
impulse response filter (FIR) design problems can be 
satisfactorily handled by using the McClellan-Parks (MP) 
algorithm [1] for weighted equiripple filters. These filters 
have the advantage of providing the designer with the most 
optimal design in the sense of smallest filter length for a 
given set of specifications. In contrast, a number of authors 
have also considered the least-squares approach for FIR 
filter design [2J-[6], As outlined in [5], [6], and [31], there 
are some advantages under certain situations where these 
methods have to be preferred over the Remez exchange 
techniques. 1 Most least-squares techniques advanced so far 
are based on solving a set of linear simultaneous equations 
by matrix inversion. 

Consider a typical low-pass design problem: we wish to 
approximate a "desired response" D(u>) with a type-1 
linear-phase FIR filter transfer function H(z) 



AT-l 



(1) 



««0 



(A type-1 filter has h n = h N _ x _ n and, moreover, the order 
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lM Remez exchange algorithm*' and the "MP approach** are used syn- 
onymously in this paper. 
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whereas the amplitude response of H(z) is [7] 

H 0 (e^) - l\^ B - W)tt = I 6„cos/ico (3) 



n-0 



n-0 



where M = (iV-l)/2 and 



The least-squares (LS) approach [2] to this problem is to 
formulate an objective function 

E^QDM-floienl— (5) 

where J? is the region 0 < o> < «r, but excluding the transi- 
tion band. Hie parameters b„ are found by rmnimiring 
E^. The actual computation of b n can be performed [2] by 
solving simultaneously a set of linear equations (or by 
matrix inversion) 

b = A' l c (6) 

where b = (b 0 b x - • • b M )'; e and A are quantities depend- 
ing upon <o p and u^. By incorporating a weighting function 
into die integrand of (5), one can attain a tradeoff between 
passband and stopband accuracies. 

Several interesting design examples can be found in [2]. 
Fig. 1 shows a typical magnitude response \H 0 (e Ja )\ of 
such a least-squares FIR filter. A particular special case of 
the above approach gives rise to the prolate-spheroidal 
wave sequence as the solution [2]. This corresponds to 
minimizing 



*Ls=f [D(*)-H 0 (e*>)] 



V 



(7) 



under a suitable constraint (such as EA* = 1). Hie resulting 
amplitude response is typically shown in Fig. 2. The pro- 
late-spheroidal solution vector h « [h Qy • • , A^-iI' (and, 
hence, b) can be found as the eigenvector of a real, 
symmetric, and positive-definite matrix [2] (which happens 
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Fig, 2. Typical magnitude response of a Kaiser window. 



to be Toeplitz). For reasonably large N, this solution h can 
be approximated to a very high degree of accuracy by 
closed-form expressions based on the zeroth-order Bessel 
function I 0 (x). Such closed-form expressions have been 
introduced and used by Kaiser [18] for the design of the 
w Kaiser window." Since the Kaiser window is a closed-form 
approximation for the above eigoi vectors, the latter need 
not be computed by elaborate eigensystem subroutines. 

Notice that the least-squares filter response shown in 
Fig. 2 is itself not a "good" low-pass response; in order to 
obtain an acceptable low-pass response (based on the LS 
approach), the passband region O^w^t), must be in- 
cluded in the objective function E^ of (5). Fig. 2 repre- 
sents a specific instance of the LS problem where the 
solution vector is an eigenvector of an appropriate real, 
symmetric, and positive-definite matrix. On the other hand, 
Fig. 1 represents another instance of the LS problem where 
the solution corresponds to an acceptably good low-pass 
response, but cannot be obtained as the eigenvector of an 
appropriate matrix. The question that arises in this context 
is, can we formulate an appropriate objective function E 
such that the filter coefficients can be obtained from an 
eigenvector of an appropriate matrix, and at the same time 
.give rise to a low-pass response as in Fig. 1? In other 
words, can we obtain a linear-algebraic generalization of 
the prolate-sequences (or the Kaiser window) so that the 
resulting vector itself has a response as in Fig. 1? 

The purpose of this paper is to address this question. 
The answer turns out to be in the affirmative, and we 
discuss the solution and several applications of this result 
Such FIR filters whose coefficients are the components of 
eigenvectors will be termed "eigenfilters" The idea of 
using an eigenvector in order to find the coefficients of a 
FIR filter has been used earlier in other contexts [16]. The 



well-known technique of Pisarenko [17] for harmonic re- 
trieval is such an example. Even though eigenvectors have 
been used in the past for filtering applications (for exam- 
ple, see [6]), we believe that the present formulation is 
novel in the sense that it takes care of the passband 
accuracy directly. Section II formulates the new quadratic 
objective function and includes design examples of such 
low-pass filters. In Section III, we show how certain time- 
domain constraints can be taken into account along with 
the usual frequency-domain specifications while designing 
these eigenfilters. Section IV presents a method for design- 
ing eigenfilters which have equiripple behavior both in the 
passband and stopband; the method is based on iteratively 
feeding bade the approximation error into the integrand of 
(5) by incorporating a weighting function. The obvious 
question that comes to mind in this context is, why should 
we use a quadratic measure at all, rather than an L m 
measure, if we are interested only in an equiripple result? 
The justification is that we can now include time-domain 
constraints if we use the eigenfilter approach, and at the 
same time obtain a more or less equiripple frequency 
response. In Section V, we show how eigenfilters can be 
used for the design of Afth-band filters (and Nyquist 
filters). The results are compared with other recently re- 
ported techniques [81 [9] for the design of Nyquist filters. 
Section VI describes how one can employ the eigenfilter 
approach to design filters which have flatness constraints 
(or tangency constraints) at any given frequency in the 
passband. In the last section, we make comments on the 
computational aspects involved in finding the appropriate 
eigenvector. 

II. Linear Phase FIR Low-Pass Eigenfilters 

Let H{z) be an FIR transfer function as in (1). We wish 
to obtain a low-pass frequency response as in Fig. 1, by 
minimizing an error measure of the form 

E = tfPv (8) 

where P 2 is a real, symmetric, and positive-definite matrix 
depending upon the design requirement, and € is a real 
vector related to h a in some simple manner (to be 
elaborated). We assume an implicit constraint th = \ to 
avoid trivial solutions. We wish the error measure E to 
reflect both the passband deviation and the stopband 
deviation from the ideal values of (2). Once such a measure 
is formed, the solution vector v is simply the eigenvector 
of P corresponding to its minimum eigenvalue in view of 
the well-known Rayleigh's principle [19]. We impose the 
additional condition that the resulting transfer function 
H(z) should have linear phase, i.e., 

where the order N-l could be either even (type 1) or odd 
(type 2) [7]. We do not consider antisymmetric impulse 

2 Notations used in the paper Bold-faced letters indicate vectors and 
matrices. Superscript / stands for transposition, whereas superscript 
dagger (t) stands Tor transposition followed by complex conjugation. 
Real symmetric positive-definite matrices P are indicated by the notation 
P*=P*>Q, and satisfy the property y'Jfy>0 for all y*Q. 

* & 
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response sequences in this section because they cannot be 
used in low-pass designs [7]. 

With H(z) satisfying (9), its frequency response takes 
the form 

H ( e J») . € -H"-WH 0 (e*) (10) 

where H 0 (e Jta ) is real-valued, given by [7] 
ht 

£ 6 w cos/iw, N-l even 

H 0 (e*) = {*J_\ . _ (11) 



( M 

Z b n cos « + » «, tf-lodd 



The quantity M in (11) is defined as M = (JV-l)/2 for 
even 7V-1, and M = N/2 for odd JV-1. Defining 



*1 ' " b M-2 



N — l even 
JV-1 odd 



(12a) 



measure for die passband can be taken as 

E , = lpe*( W )rf W = i/V(l-*)(l-r)W« (18) 

which can be written in the form 

E t = Vt p b (W) 

where P p is given by 

P=-r(l-c)(i-c)'d a (20) 

and is a real, symmetric, positive-definite matrix (unless 
» = 0). Thus, the total measure to be minimized is 

E = VPb (21) 

where 

P=(l-a)P p + aP s , (22) 
The quantity a, which is in the range 0 < o<l, controls 



and 



*(«) 




cos a cos2w : -cosAfw]', N — 1 even 

|M-iJ<o]', JV-lodd 



3a> 



(12b) 



cos 



-cos 



the relative accuracies of approximation in the pass and 
stopbands. Notice that the dements of P are given by 

(1-a) f» 

P(n,m) = / (1 - cos nw)(l - cos mu) 



a 



we can write (11) as 

(13) 

For notational simplicity, c(«) will often be denoted as r. 
For even N - 1, the quantities are as in (4). For odd 
N-\ 9 b n = 2h M _ x _ tt . 

With the "desired response" as in (2), the "stopband 
error" can now be defined as 

E s ~- r[D{u)-H o {e*)] 2 d0-^ C Vcc'bdu^VPsb 

V J *>S m **S 

(14) 

where P s is given by 

P s = - fee' da (15) 

IT J us 

and is a real, symmetric, and positive definite matrix 
(unless to 5 = tt, which is a case of no interest). 

If the passband error measure E p is also defined accord- 
ing to the integral of (5), the total error measure cannot be 
written in the form (8). Accordingly, let us define E p 

differently. First, notice that the zero-frequency response is eigenproblem. Given the band edges co p and u> S9 and the 
given by 

H 0 (e'°) = Vb (16) 
where the vector 1 is defined as 

1 = [11 ••.!]'. (17) 



+ — f r (cosnoicosmu>)dio 9 Nodi (23a) 

p( B , w) =^jr(i-cos(»4) tt ) 

cos|/n + — jwj du> 

+ ^h( n+ i) w ]h( m+ i)T w ' 

AT even. (23b) 

In summary, we have been able to formulate the linear- 
phase low-pass FIR design problem in the form of an 



The quantity e p {o) = (l-cyb 9 therefore, represents the 
deviation of the response H 0 (e J °) from the zero-frequency 



parameter a, the matrix P can be computed* It is easy to 
obtain closed-form expressions for the integrals in (23), 
and, hence, the elements P(n,m) are easily computed 
once a,«p and w s are known. It then remains only to 
compute the eigenvector of a real, symmetric, and 
positive-definite matrix corresponding to the smallest 
eigenvalue. The resulting filter is guaranteed to have linear 



response. 3 Accordingly, a positive-valued (quadratic) error phase because the vector b rather than the vector h is 

directly involved in the optimization problem. The eigen- 
3 The only motivation for taking zero-frequency as a reference for vector b can be used to obtain the filter coefficients h n in 
passband error formulation is that it brings the vector b into the ry\ \y e QO W proceed with design examples to demonstrate 
reference, and this enables us to write £_ as a quadratic in 6. This in turn v '* r 
leads to the eigenformulatton. the procedures. 
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Fig. 3. (Ex. 1) Magnitude response plot of Kaiser window using eigen- 
filter and /(/*) approximation. 
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Fig. 4. (Ex. 2) JJa) Magnitude _resgons£j>lot of FIR eigenfilter. 



details for the FIR eigenfilter. 



Example 1: Let « s = 0.05«r and let a=l in (22). Then 
to does not enter into the optimization problem. With 
AT - 1 = 30, the eigenvector b corresponding to the mini- 
mum eigenvalue of P was computed The resulting 
frequency response H 0 (e J ") is shown in Fig.. 3 in solid 
lines. The broken line represents the transform of the 



<0, + CO, * 7T 




Fig. 5. Typical amplitude response of a half-band FIR filter. 

Kaiser window for a window length of 31, and an ap- 
propriate value of the window parameter [7] given by 
P = 2.0. This example is only a .demonstration that the 
Kaiser window is an excellent closed-form approximation 
for certain eigenfilters (the prolate spheroidal sequences). 

Example 2: A low-pass filter with #-1 = 28, w p = 
0.3*r, w s = 0.4*7, and a = 0.1 was designed using the above 
approach. The resulting frequency response is shown in 
Fig. 4, which also includes the plot for the case of a = 0.5. 
The effect of a is clearly seen from the figure. 

Comments on the Choice of a: It is clear from (22) that a 
larger value of a leads to better stopband attenuation at 
the cost of increased passband ripple. Given the set of 
specifications <o p , co 5 , passband tolerance S x and stopband 
tolerance S 2 , it is necessary to estimate N-l and a in 
order to design the eigenfilter. An approximate estimate 
for N-l can be obtained based on the relations in [20]. 
Even though the estimates in [20] hold only for equiripple 
designs, the required order for an eigenfilter is only slightly 
larger. The choice of a governs the ratio o^/Sj. At this 
time, we do not have a procedure for estimating a, starting 
from a desired 8i/8 2 ratio - Further study of the behavior 
of eigenfilters is necessary in order to fill this gap. 

The following example demonstrates the design of a 
half-band filter [8], [11] based on the eigenfilter approach. 
For this example, the quantity a does not come into the 
picture at all, in view of the indirect procedure described 
below. 

Example 3: A linear-phase FIR half-band filter is a 
type-1 low-pass filter having an amplitude response 
H 0 (e J ") that is symmetric with respect to v/2 as shown in 
Fig. 5. Note that + = tt, and 8 X = fi 2 f° r such filters. 
It can be shown that for such filters h n satisfies 



0 



for n- 



N-l 



= even * 0. 



(24) 



A direct design of such filters based on any design meth- 
odology (such as the McClellan-Parks program, or eigen- 
filter methods) invariably leads to a filter which does not 
satisfy (24) exactly, due to the computational inaccuracies. 
In addition, since about half of the impulse response 
coefficients are known anyway, it is only judicious to 
eliminate the known variables from the design problem so 
as to reduce computational time. This can be accomplished 
by designing H(z) tjy a two-stage process. The first step is 
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Fig. 6. Amplitude response of G(z). 
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Fig. 8. (Ex. 4) (a) Typical magnitude response of a bandpass FIR 
filter, (b) Magnitude response plot of a bandpass FIR eigenfflter. 

the specifications 2u> p and B in Fig. 6 are known. The 
order of G(z) is half the order of H(z) and G(z) can be 
designed much more quickly. Moreover, the computation 
of H(z) as in (25) ensures that (24) is satisfied exactly, 
with no error. More details in this connection can be found 
in [32]. 

The eigenfilter approach can be used easily to design the 
odd-order one-band filter G(z) with response as in Fig. 6, 
simply by taking a = 0 in (22) and letting 



IORRRLIZED FREQUENCY 

(*>) 

Fig. 7. (Ex. 3) (a) Response of the half-band ei&enfHter. (b) Passband 
details for the half -band eigenfilter. 



the design of a one-band linear phase filter G(z) of odd 
order(W - 1)/2 with response as in Fig. 6, where 8 = 25. 

Hie passband of G(z) is the interval 0 < o> < 2u> p9 and 
the region 2o> p < v is the transition band. Since (N— 
l)/2 is odd, G(e J ") is automatically zero at <o = *r [7]. 
Now, H{z\ which has a desired response as in Fig. 5, can 
be obtained as . 

G(z 2 ) + z-<»-W , x 

H(z)= 1 2 . (25) 

Given the specifications <o p and 5 for the half-band filter, 



,-*-I/*v<l- e Xl-«)' 



die. 



(26) 



The eigenvector b contains information about the impulse 
response g„ of G(z). The impulse response of H(z) is 
found as 



0.5g n/2 , n even 



n odd # 



(27) 



1 



AT-1 



Fig. 7 shows the frequency response of a half-band filter 
designed using the above eigenfilter approach. Here, N — l 



« 34 and <o p = 0.4225*r. 
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Example 4: The design of high-pass and bandpass 
eigenfilters can be accomplished with equal ease. For 
example, a type-1 bandpass filter can be designed by 
defining the objective function E as E = aE l + PE 2 + (1 — 
a-fi)E 3 , where a,p>0 and l-a-/*>0. E x and E 2 
represent the stopband errors 

E^-b' Pec' dub E 2 = -vC ce'dub 

and £ 3 is the quadratic measure of passband error defined 

£ 3 - »*' P(a - c)(a - c) ' die ft. 

Here, a is a constant vector of the form 

a = [1 cos<o 0 cos2<i) 0 - — J 

and co 1? <o 2 , co 3 , and « 4 are the band edges as shown in 
Fig. 8(a). The quantity <*> 0 is taken as 6J 0 = (o3 2 + co 3 )/2. 
Fig. 8(b) shows the response for a design example where 
the order is AT - 1 = 50, and the band edges are = 0.3tt, 
« 2 — 0.35*r, « 3 = 0.7sr, and « 4 — 0.8tt. 

iii. elgenfilters with simultaneous 
Time-Domain Constraints 

Error functions that are more general than (21) can be 
formulated in order to take into account time-domain 
constraints. As an example, let us assume that we wish to 
design a linear-phase low-pass FIR filter as in Fig. 9. Let 
us assume that, in the input signal x(n\ there are occa- 
sional segments of an unwanted waveform, i.e., a finite 
duration waveform of the type 

s 0$ s l9 "^s L _ v (28) 

The time of occurrence of this L-point waveform is un- 
known, but its shape is known. We would like H(z) in (1) 
to be such that, in addition to being a good low-pass filter, 
it minimizes the effect of s k on the output signal y n . Thus, 
let 



«0 f l,---,L + JV-2 (29) 



represent the (L + N — l)-point sequence at the output in 
response to the L-point input in (28). In matrix form, we 
have 



where - 
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(30) 



*<n>+s(n) 




y(n)+r(n) 



Fig. 9. System with input noise. 



r=[r 0 ^ • ^ +J ,_ 2 ]' (31) 
A = [A 0 h 1 -h„_ x )'. (32) 

A quadratic measure of the output error is given by r'r. 
Define a normalized measure of output "noise" by 

h'S'Sh 



E N = 



(33) 



where 



In order to preserve the linearity of the phase of H(z% it 
is necessary to write E N in terms of the b vector so that 
the b vector can be made the "unknown" vector of the 
optimization problem. For this, notice that the vectors h 
(eq. (32)) and b (eq. (12a)) are related as h = Db, where 



1 



and 



(0 0 0 

6 6i 

0 10 

2 0 0 

0 10 

0 0 1 

[0 0 0 

'0 0 0 

0 0 0 



1 

"-2 



0 
0 
0 
0 
0 

il 



0 1 

1 0 



for AT -1 even (35) 




0 0 0 
0 0 0 



Thus, E N simplifies to 
where 



for AT -1 odd. 



1 0 
0 1) 



E N = VP N b 



D'S'SD 



(36) 
(37) 

(38) 



The overall quadratic objective function to be minimized is 
therefore 

E = (l-a-P)E p + aE s + pE N 

- A'[(l - a- p)P p + aP s + pP N ]b (39) 
where a > 0, ft > 0, and a + P < 1. Once again, for a given 
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Fig. 11. (Ex. 5) (a) Frequency response of an eigenfilter having simulta- 
neous Ume-domain constraint (b) Passband details. 

set of parameters (a p9 o> s , a, and /?), the filter design can 
be completed by finding the eigenvector of the matrix 

P=(l-a-fi)P p + aP s + pP N (40) 

corresponding to its minimum eigenvalue. Notice that 
P = P'>0. 

Example 5: As a numerical example, let #-1 = 28, 
L - 1 = 19, = 0.3(*r), a> s = 0.4(ir), and a = 0.5, and let 
the waveform s n be as in Fig. 10. Fig. 11 shows the 
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Fig. 12. A typical step response for a low-pass filter. 

resulting frequency response of the linear-phase eigenfilter 
with 0 = 5X1O" 4 (solid lines) and with 0 = 5.OxlO~ 3 
(broken lines). 

In an analogous manner, it is possible to put constraints 
on the unit step response. Notice that the first k+l 
coefficients of the step response t n are given by 
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(41) 
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It is often required to control (minimize) the energy in the 
step response over the interval 0 < n < k y where k is the 
time taken for the beginning of a monotone rise (Fig. 12). 
Once again, the energy t't can be written in the form 
VPjb, where P T = Ff > 0. Thus, additional quadratic con- 
straints can be introduced into the problem in this manner 
before computing the eigenfilter. 

IV. Equiripple Linear-Phase Eigenfilters 

There exists a powerful technique [1] based on a Remez 
exchange procedure for designing equiripple filters by 
minimizing the norm of the response error. In view of 
the availability of such excellent techniques, the question 
arises as to why we should manipulate an eigenfilter to 
give an equiripple response. The answer lies in the fact that 
eigenfilters can be designed so as to take into account 
other design constraints such as time-domain constraints, 
as illustrated in Section in. In addition, eigenfilters can be 
used for designing Nyquist filters by elegant modifications 
of the positive-definite matrices appearing in the objective 
function (Section V). For such applications, if we have a 
technique for equalizing the error over each band, the 
result is of a more or less equiripple nature, and the peak 
error tends to get minimized. With this motivation, we 
describe the algorithm for accomplishing this as follows. 

The quantities E p and E s in (18) and (14) are the error 
measures in the frequency domain. Suppose we incorpo- 
rate a positive-valued weighting function W(w) into the 
integrands 

Es = ~f el(»)W(»)d* (42) 
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and 



where 



E p = ^£*e$(<*)W(»)d<* (43) 



e s («)-tt *,(6>) = (l-c)V 



(44) 



For uniform weighting (Le., with W(w) = l), the errors 
eJto) and e s (o>) of the resulting eigenfilter are large near 
the band edges and tend to fall off at points away from the 
band edges (see Fig. 4, Example 2). Accordingly, if W(u>) 
is larger near the band edges than at other points, the error 
tends to get equalized. Since the exact shape of the eigen- 
filter's error responses e p (w% e $ (u) is not known a priori, 
we follow an iterative procedure to identify the ap- 
propriate W{w) as follows. 

Let * p ,jfc(«) and e s%k (<*) be the error curves for the 
eigenfilter at die end of the &th iteration. Then the weight- 
ing function for the (k + l)th iteration is taken to be 



f^*(«H*,.*(«)l. 



0 < « < «p 
6>5 < w < m 



(45) 



Notice that either the errors e pk (<*>) and e 5aJb (a) or their 
envelopes (which are monotone) can be used in (45). Once 
the errors become equiripple, the weighting function does 
not alter any more with further, iterations, and we can stop 
the iteration. Note that, once we employ such weighting 
functions in (42) and (43), the integration should be per- 
formed numerically, thus increasing the design time. 4 
Another possible scheme for finding would be to 

replace e, v *(oX*s,*( tf ) m ( 45 ) a low-order poly- 
nomial that fits the envelope of the error. This would avoid 
the need for numerical integration, but should be weighed 
against the time required to obtain the polynomial fit 

Example 6: The eigenfilter considered in Example 2 
was redesigned using the above procedure. The resulting 
design after 16 iterations has frequency response as shown 
in Fig. 13 (solid lines). The figure also shows the response 
of an optimal filter designed using the Remez exchange 
algorithm [1] (broken lines) with the same a> pJ o> s and the 
same order N-l as the eigenfilter. Notice that, for the 
equiripple eigenfilter, S p = 0.032 and S 5 = -31.25 dB, 
whereas, for the Remez-exchange-based optimal filter, S p 
= 0.031 and S s = -31.2 dB. Indeed, the tolerances 
achieved by the eigenfilter are very close to optimal. The 
numerical integration involved in the updating equation 
(45) was performed by using Simpson's rule, with 100 grid 
points. As a rule of thumb, the number of grid points can 
be taken to be about 20 to 30 times the approximate 
number of extrema in the frequency band under considera- 
tion. We are not aware of a proof that the above technique 
for equiripple designs always converges to the optimal 
solution. However, all of our design examples demonstrate 
such convergence. 

4 We learned from a reviewer that this kind of error feedback technique 
is in fact known [27], [28]. 
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Fig. 13. (Ex. 6) (a) Magnitude response of equiripple FIR filter de- 
signed using eigenfilter and Remez exchange algorithms, (b) Passband 
details. 



V. Design of Nyquist Eigenfilters 

The design of digital Nyquist filters for generating a 
band-limited pulse for data transmission with minimum 
intersymbol interference has a significant role in the design 
of digital modem systems [12]-{15]. The impulse response 
of a Nyquist filter of length N satisfies 
/ N-l\ 

/i n = 0 for In — I = nonzero multiple of K. 

These have also been known as ATth-band filters [8]. Such 
filters are also of great interest in the design of interpola- 
tion filters [31, ch. 4] in multirate signal processing. The 
cutoff frequencies co p and <o s of a Nyquist filter satisfy the 
condition w p + « s = 2flr/^ where AT is the intersymbol 
duration. 

The design of Nyquist filters requires an algorithm that 
incorporates both time-domain and frequency^domain 
constraints. Such designs using linear programming or 
based on the Remez exchange algorithm have recently 
been suggested by several authors [8], [9] T [IS]. In this 
section, we would like to apply the eigenfilter approach, 
with some additional modifications on the objective func- 
tion, to design the Nyquist filter. Recall from Section II 
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Fig. 14. (Ex. 7) (a) Nyquist dgenfflter response, (b) Passband details. 

that the error E is given by (21) for a low-pass design, 
where P is as in (22). We note that every if th coefficient 
except the mid-term of the impulse response of a Nyquist 
filter is exactly zero; hence 

*V = [b 09 b lt • • • , 6jr-i»0, b K + u • • • , b 2K _ l9 0, b 2K+u • • • ] . 

(46) 

Therefore, those rows and columns of P whose indices are 
multiples of K do not contribute to the total error E = tfPb. 
To take advantage of this fact, we restate our problem as 
minimizing the error J?j * tVF'tY, where 
*" = [b Qf b l9 - • • , b K - l9 * * , b 2K -i 9 *2ic+i» ' 
and 



Note that the order of P f is smaller than the order of P; 
we, therefore, save time in the design process. Once we 
find the appropriate eigenvector of P\ we insert the zeros 
that are required in the Nyquist filter impulse response. 
This completes the design. Once again, we can feedback 
the error in both passband and stopband of the previous 
iteration to compute P and P' according to (42) and (43) 
of Section IV to obtain an approximately equiripple 
frequency response for the Nyquist filter. (In such exam- 
ples where we have both time- and frequency-domain 
specifications, it is, in general, not possible to get an 
exactly equiripple solution.) Notice that the resulting filter 
automatically has linear phase because we are directly 
solving for the vector b. 

Example 7: A Nyquist filter with order #-1 = 38, a p 
= 0.15(«r), « s = 025(flr), JSC = 5, and a = 0.95 is designed 
using the eigenfilter approach. Hie frequency response of 
the Nyquist filter (solid lines) as well as of the equiripple 
Nyquist eigenfilter (broken lines) obtained from the tech- 
niques described in Section IV are shown in Fig. 14. Note 
that the transition band of the equiripple Nyquist eigen- 
filter is smaller than that of the Nyquist eigenfilter, and 
explains why the peak errors are larger for the equiripple 
design. 

Example 8: A FIR linear-phase Nyquist eigenfilter with 
length 39, c0p = 0J125(sr), % = 0.2875(*r), a = 0.98, and 
K = 4 is designed. This is compared to the Nyquist filter 
with the same specifications, designed by linear program- 
ming [91 and also the design based on the Remez exchange 
algorithm [8]. We summarize the passband error 8 p and 
stopband error 8 S for these approaches. 



] (47) 



Eigenfilter Remez 
8 p (db) -0.45 -0.335 
8 s (db) -3321 -29.7 



Linear Programming 
-0.45 
-33 



Hie plot of the frequency response using the eigenfilter 
(solid curve) and MP approaches [8] (broken lines) are 
shown in Fig. 15. The frequency-response plot of the 
Nyquist filter designed using linear programming can be 
found in [9]. Notice that the eigenfilter response is slightly 
better in the stopband than that of the Remez exchange 
method, and slightly worse in the passband. A readjust- 
ment of a in (22) can produce results much closer to the 
Remez exchange procedure. 



r o,o 



Or-1,0 



r 2jr-i f o 



M),JC-1 



M>,X>1 

P i.jr+i 



^2jc+i.#:+i 



MUX-l 



? 1,2JT+1 



^21C+l,2Jr-l 



2JT+1 



r JC+l,2*+l 



^2JC-1 



,2X*+1 



r 2JT+l,2Jt+l 



(48) 
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Fig. 15. (Ex. S) (a) Magnitude response of Nyquist filter, designed 
using the eiacnnlter approach and Remez exchange algorithm, 
(b) Passband details. 



and 



this section, we impose the flatness constraints at any 
frequency in the passband using the eigenfilter approach. 

Consider an even-order filter (type 1) H{z) of order 
AT — 1. Its amplitude response can be written as a function 
of sin 2 («/2), Le., 



H 0 (e J ") = E *„cos««= £ c„cos 2 -(^) 

„«o - «-0 



Af 



(49) 



Suppose the first L +1 coefficients (except rf 0 ) of H Q (e J ") 
are zero, i.e., 



+ d u sin 2M 



(!)• 

(50) 



Clearly, the first 2L + 1 derivatives of H 0 (e Jia ) are zero at 
« = 0, Le M 



-Ti;H 0 {e*>) = 0, 1 < k < 2L + 1 . (51) 



This fact suggests a procedure to design filters with flatness 
constraints using the eigenfilter approach. Defining 

d=[d 0 d x d 2 - d M . x d M V N -I even (52) 



sin 



,2(Af- 



we then have 

JZo(e>") 

To achieve a 2L + 1 degree of flatness at <o = 0, we define new sequences J' and t>' as 

<f'=»[<f 0 d L+l d L + 2 ~' d M-i d i*V 



t/=[l sin^ +1 >(y) sin^ +J >(y) sin*"-^^ sin 2 "(y)]'. 



(53) 

(54) 

(55) 
(56) 



VI. Eigknfilters that have Arbitrary Degree Incorporating this ^C^aifc™ fcenew 
of Flatness at Any Given Frequency in the objective function E = (d')'Fd', where i> is as in (22) with 



Passband 



Digital FIR filters with a maximally flat frequency 
response at frequency co = 0 have been studied by 
Herrmann [21], Kaiser [22], and Vaidyanathah [23]. FIR 

filters of this class find applications in several filtering n 0 

problems, as pointed out by Kaiser and Stdglitz [24]. In where 1' = [1 0 .0 



P s ~-C Jifdu 

p p =-r(r-*')(r-*'Yd<* 



(57) 



(58) 

0]'. Now if we find the eigen- 
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vector corresponding to the minimum eigenvalue of P, the 
resulting response has the desired flatness at a> = 0. We 
can easily generalize this approach when the flatness con- 
straints are required at an arbitrary frequency. Thus, the 
frequency g> 0 at which the first 2L + 1 derivatives are zero 
can be varied by expanding H Q (e JtJ ) in terms of sin 2 " ( w - 
« 0 /2)in(49). 

VII. Eigenvector Computation and Related 
Issues 

A major part of the design time for our method is spent 
in the computation of the eigenvector. Since we are inter- 
ested in only one eigenvector (corresponding to an ex- 
tremal eigenvalue), this computation can be done effi- 
ciently (without invoking general methods such as the QR 
technique [29]). It is well known [25], [26] that to compute 
the dominant eigenvalue and its corresponding eigenvec- 
tor, the iterative power method is simple and fast if the 
ratio l^tf/^Ar-il is large where \ k are the eigenvalues of P 
with 

l*»l>l*»-d>"- >l*al>l*d- 

At the k + 1th iteration of the power method, a vector 
x k+l is computed from the previous iterate x k as 

yk+x = P*k (59a) 

** + i-jWIIJfc + ill (59b) 
where denotes the L 2 norm of vector v. The difference 
between x k and x k+l , defined by \\x k + x — x k \\ is compared 
to a prescribed small constant c, and if 

(60) 

then x k+l is a good approximation of the eigenvector 
corresponding to the maximum eigenvalue. A typical value 
for c is about 1.0 X 10~ 6 . 

We, however, wish to compute the minimum eigenvalue 
and its corresponding eigenvector. In other words, at each 
iteration, we would like to compute 



(61) 



Given x kJ we would like to find x k+l without inverting 
P~\ Rewrite (61) as 

** = *** + i (62) 

It is well known [25], [26] that a real, symmetric, and 
positive-definite matrix P can be decomposed into 

P = LV (63) 

where L is a real lower triangle matrix. Equation (62) then 
becomes 



Let 
Then 



x k -LL'x M . 

X k= Lv k+V 



(64) 
(65) 
(66) 



We can find v k+1 in (66), given x k and £, by recursively 
solving a set of linear equations. Let be the element in 
the ith row and yth column of L; then we can show that 

"* + i00 = t-(**(*) - "LI.jO M U))- (67) 

Since P is positive-definite, l nn in (67) are evidently non- 
zero. Using (67), we first solve for o ft+1 (0) and solve 
recursively for all v k + x (n) for l^n^N — 1, where N is 
the dimension of L. (In (67), x k (n) denotes the nth 
element of vector x*.) It takes N(N - 1)/2 multiplications, 
N divisions, and N(N-\)/2 additions to compute v k + v 
Similarly, we can find x k+l in (65) given v k ^ l and L. Hie 
total time required per iteration is thus 



('« + '*) 



N(N-l) 



X2 



where t a> t m > t d are, respectively, the required computer 
times for addition, multiplication, and division. 

Note that the speed of convergence depends on the ratio 
X 2 /Ai since we are dealing with P~ l rather than P. 
Instead of proceeding as in (65) and (66), one could invert 
P beforehand, store it as Q, and then perform the iteration 

y k+l = Qx k (68a) 

**+i-jWIIJfc + ill- (68b) 
The operation (68a) requires N z multiplications and N(N 
— 1) additions, and the time required for this is 

f fl #(JV-l)+*JV 2 

which is nearly the same as the time required for per- 
forming (65) and (66). However, there is an overhead cost 
associated with the computation of Q. 

Example 9: Table I indicates a comparison of design 
time for half-band filters, using the eigenfilter approach 
and the Remez exchange approach, for various tolerances. 
The ratio X 2 /X x the number of iterations required for 
eigenvector computation are also given. We observe that 
X 2 /Xi is large and, hence, the number of iterations for 
eigenvector computation is impressively small for all the 
entries in Table I. 

We also summarize below the design time, number of 
iterations, X 2 /\ 1 , 8 V and S 2 for the Nyquist filter de- 
signed in Example 8 (solid curve in Fig. 15): 



CPU time Number of X 2 /Xj 
(sec) iterations 
3,9 3 792 



«! (db) S 2 (db) 



-0.45 -33 



Our experience based on a large number of design exam- 
ples has convinced us that the ratio Xj/Xj is large in all 
practical cases. Accordingly, the eigenvector computation 
never creates any numerical or stability problems and is 
invariably fast Single precision arithmetic is found to be 
sufficient in all design examples. 
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TABLE I 

Comparison of the Eigenfilter and the Remez Exchange 
Approaches in the Half-Band Filter Designs 



eigenfilter 



REMEZ EXCHANGE 
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j 1 Norn! 
item 


bcr of CPU time 


If-1 
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CPU 
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4 0.8 
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.425* 


tt 


0u0317 


309 


3 0l5 


IB 


0l0365 
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n-l is the required order for peak 
width 2(tr/2- a p ). 



passband ripple &i and transition band- 



There exist some recent methods for computing eigen- 
vectors (corresponding to extremal eigenvalues) based on 
gradient techniques [30]. These could prove to be even 
faster than the iterative power method, but we have not 
studied this possibility in the context of our paper. 

vm. Concluding Remarks 

In this paper, we have formulated linear-phase FIR 
design problems in the form of eigenproblems. The filter 
coefficients are computed from eigenvectors of certain 
matrices which represent the specification requirements. A 
set of specification requirements directly leads to the 
Kaiser-window frequency response. Eigenfilters can be 
designed in such a way that various time-domain require- 
ments are simultaneously taken into account By means of 
an error-feedback mechanism, an iterative method also has 
been presented that leads to near-equiripple eigenfilters. 
The eigenfilter approach is also used to design Nyquist 
filters and filters with flatness constraints in the frequency 
domain. In conclusion, a wide variety of specifications can 
be taken into account by the eigenfilter approach. Hie 
design examples presented and the design times compare 
favorably with several other methods reported in the litera- 
ture. 
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